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We discuss the lattice formulation of gauge theories with fermions in arbitrary representations of 
the color group, and present in detail the implementation of the HMC/RHMC algorithm for simulat- 
ing dynamical fermions. We discuss the validation of the implementation through an extensive set of 
tests, and the stability of simulations by monitoring the distribution of the lowest eigenvalue of the 
Wilson-Dirac operator. Working with two flavors of Wilson fermions in the adjoint representation, 
benchmark results for realistic lattice simulations are presented. Runs are performed on different 
lattice sizes ranging from 4^ x 8 to 24^ x 64 sites. For the two smallest lattices we also report the 
measured values of benchmark mesonic observables. These results can be used as a baseline for 
rapid cross-checks of simulations in higher representations. The results presented here are the first 
steps towards more extensive investigations with controlled systematic errors, aiming at a detailed 
understanding of the phase structure of these theories, and of their viability as candidates for strong 
dynamics beyond the Standard model. 



I. INTRODUCTION 

Recent algorithmic progress in simulations of QCD 
with dynamical fermions have shown that current com- 
puting facilities can reach beyond the quenched approxi- 
mation without compromising the robustness of the for- 
mulation [ll-Q. Besides their great relevance for lattice 
QCD, these improvements allow numerical simulations 
to be used as a quantitative, nonperturbative tool to 
study more generic strongly-interacting theories where 
the role of the fermion determinant is crucial: exten- 
sive Monte Carlo studies will enable to tame the system- 
atic errors and therefore to obtain reliable information 
about the strong dynamics of theories beyond QCD. Ex- 
amples that have a direct impact on physics beyond the 
Standard Model (BSM) include technicolor models 0, |i| 
(see d, for reviews), orientifold theories [lH, and su- 
persymmetry (see [l^ for a recent review on the nonper- 
turbative formulation of supersymmetry) . Moreover the 
theoretical progresses in formulating QCD on the lattice 
provide the necessary tools for rigorous studies of BSM 
nonperturbative phenomena. This is of paramount im- 
portance in order to obtain reliable results as we move 
into the unknown territories of strongly-interacting BSM 
theories: the clean field-theoretical framework developed 
for QCD needs to be exported to the studies of BSM 
strong dynamics in order to obtain lattice results that 
are trustable and therefore can have an impact on the 
LHC phenomenology. 

Simulations of the theories listed above share a com- 
mon technical challenge: they all require fermion fields 
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in representations larger than the fundamental, gauge 
groups SU(A'^) with generic number of colors N, or a 
combination of both. Therefore a central ingredient of 
these simulations is a generalized Dirac operator, that 
can act on vectors of arbitrary dimension in color space. 
We developed such an operator and tested its inversion 
in a recent study of the quenched meson s pec trum in the 
large- limit of pure Yang-Mills theories [i3[ . The same 
computational framework has also been used to assess 
the validity of the planar orientifold equivalence by nu- 
merical simulations in the quenched approximation fl^. 

Theories with two Dirac fermions in the two-index 
symmetric representation and small number of colors 
{N — 2,3) are currently under scrutiny as viable can- 
didates for strongly-interacting BSM candidates [T5l - [l7| . 
An important feature of these theories is that the num- 
ber of fermions is close to the critical value where an IR 
conformal fixed point may appear [l8| . The vicinity of 
the IR fixed point can induce a walking behaviour of the 
running coupling, which would make such theories phe- 
nomenologically viable. Whether or not this is really the 
case beyond perturbation theory can only be established 
by numerical simulations of such theories defined on the 
lattice. Recent studies of the running of the renormalized 
coupling in the Schrodinger functional (SF) scheme [l^ 
with fundamental staggered fermions have provided first 
evidence of the existence of a conformal fixed point for 
Uf = 12 A similar study has also appeared for the 
SU(3) theory with fermions in the two-index symmet- 
ric representation [2l|, and for the SU(2) theory with 
fermions in the adjoint representation [2^, [23| . The re- 
sults presented in those works confirm that the SF is 
indeed a valuable tool to study the RG structure of the 
theories under consideration. Further studies on the ex- 
istence of an IR fixed point with fundamental fermions 
have appeared in the literature p4| . Several complemen- 
tary approaches will be needed in order to fully under- 
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stand numerical simulations in the vicinity of an IR fixed 
point. 

In this paper we discuss in detail the implementation of 
an algorithm for simulating gauge theories with dynam- 
ical Wilson fermions in an arbitrary representation, and 
arbitrary gauge group S\]{N), and present the results of 
preliminary runs which are useful as a benchmark for fu- 
ture simulations. At this stage in our investigations, the 
emphasis is on the algorithm rather than the optimiza- 
tion or the phenomenological results. Sect. HI] summa- 
rizes the notation used for the Dirac operator, and intro- 
duces the conventions to deal with the higher representa- 
tions. The HMC algorithm for two flavors is described in 
Sect. IIII Al In order to be able to simulate theories with 
an arbitrary number of flavors, the RHMC algorithm Q 
has also been developed, as described in Sect. IIIIBI Let 
us emphasize once again that, as numerical simulations 
move away from the well-known realms of QCD, we find 
it necessary to have a robust field-theoretical formulation 
and a tight control over the behaviour of the numerical 
simulations. These are necessary conditions if numerical 
simulations want to be in a position to test candidate 
theories for strong electroweak symmetry breaking with 
any degree of confidence. 

We study the behaviour of the algorithm at difi^erent 
points in the space of bare lattice couplings by moni- 
toring the time-history, the probability distribution and 
the correlation time of the plaquette and the lowest eigen- 
value of the Dirac operator, in Sect. II VI In the same Sec- 
tion we also study how the force is split in the terms of the 
rational approximation, in analogy with Ref. j0| where 
fundamental fermions were considered. We present re- 
sults for 43 X 8, 8^ X 16, 12^ x 24, 16^ x 32, 24^ x 64 lattices, 
for the SU(2) theory with two flavours in the adjoint rep- 
resentation at two values of the coupling constant /? = 2 
and P = 2.25. The emphasis in this work is on the details 
of the implementation, and on the behaviour of the algo- 
rithm for the theory under consideration. These results 
are intended as benchmark for other softwares that sim- 
ulate the same theory. When the exactness of numerical 
simulations cannot be judged by comparing with experi- 
mental data, a cross-check of all the available simulation 
softwares is essential before debating about the interpre- 
tation of the results. Up to date this preliminary cross- 
check is still missing, although several more complex and 
phenomenologically relevant results are currently avail- 
able in the literature f2l[2l[25l430t (for studies with dif- 
ferent number of colors or fermionic representation, see 
IM 0, HI mi-Iii] ) . Therefore the authors encourage dif- 
ferent groups to publish analogous studies, in order to fill 
this gap. 

Some phenomenological results for the simulations pre- 
sented in this work have been already published in a short 
paper [sOj (see also [5l|, [S^ ) ■ A more detailed discussion 
of these results is currently in preparation. 



II. WILSON FERMIONS IN A GENERIC 
REPRESENTATION 

In four-dimensional Euclidean space the massless 
Wilson-Dirac operator is written following the notation 
in Ref. [H: 

where V and V* indicate respectively the forward and 
backward covariant derivatives: 

V^/(x) = i pix, fi)fix + /i) - fix)] , (2) 

v;f{x) = ^[f{x)~uix-^i,^lyfix-^,)] , (3) 

and U (x, fi) denote the link variables as usual. This ex- 
pression is easily generalized to an arbitrary representa- 
tion R of the gauge group; the action of the massive Dirac 
operator on a spinor field ip yields: 

D„^^{x) = {D + mo)ij{x) = 
= (4/a + mo)'4'{x) — 

+ {l + J^)u''{x-^l,^ly^{x-fl)}, (4) 

where are the link variables in the representation R, 
and mo is the bare mass. 

Let Tp, be respectively the generators in the funda- 
mental representation and in a generic representation R. 
Our conventions for the normalization of generators and 
other group-theoretical factors can be found in App. El 
The explicit form of the generators used in this work is 
given in App. |B] If the link variables in the fundamental 
representation are written as: 

U{x,fi)^exp[iuj''{x,n)T^], (5) 

then the link variables in the representation R are given 
by: 

U''{x,fi)^exp[zu;-{x,fi)T^], (6) 

where the functions a;''(a:;, /i) are the same in both equa- 
tions. The explicit relation between link variables in dif- 
ferent representations can be worked out explicitly for 
each individual case (see App. |B]). For instance for the 
adjoint representation of SU(iV) we have the well-known 
formula: 

K'^ = ^tr [T^UTl^U^ . (7) 
The fermionic lattice action: 

Sf = a''Y,^i^)D,ni^i^) (8) 

X 
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is quadratic in the fermionic fields and therefore its con- 
tribution to the partition function can be computed ex- 
actly, yielding the fermionic determinant, det Dm- In 
numerical simulations, the latter determinant is conve- 
niently replaced by the determinant of the Hermitean 
Dirac operator Qm, defined as usual as: 

Qm = 75-D,„. (9) 

The fermionic determinant for an even number ri/ of de- 
generate flavours can be represented by introducing com- 
plex pseudofermionic (bosonic) fields (f>,(f>^: 

(det - (detQ,„)"' = J 'D^'D<f>^ e'^f , (10) 

where the pseudofermionic action has the generic form: 
5/ = a4^0t(^)(Q2^)-'0(3,)^ l = nf/2. (11) 

Note that we have used the same symbol 5*/ to denote the 
fermionic action written in term of the usual Grassmann 
variables and the pseudofermionic action written in terms 
of bosonic fields that is used in the numerical simulations; 
the meaning of the symbol will be specified by the context 
in which it is used. In the above expression the square 
of Qm is used and not the matrix itself because Qm may 
not be positive definite. The sum over spin and color 
indices in the above formulas has been omitted in order 
to simplify the notation. 

For the gauge action we always use the Wilson action 
written in terms of elementary 1x1 plaquettes built of 
link variables in the fundamental representation: 

=/^E (l-^I^^^t'-^M-) , (12) 

fj.<i' ^ ^ 

where f3 = 2N/ is the lattice bare coupling. The algo- 
rithm can be generalized to other gauge actions without 
any major difficulty. 

The partition function for the theory with n/ flavors 
of dynamical quarks can therefore be written as: 

Z = j 2?f/2?,/)2?0te-^«(^)-^/(^.^.^'). (13) 

It is also possible to work exclusively with link variables 
in the higher representation i?, by writing the gauge ac- 
tion as a function of the latter. We expect the two choices 
to yield the same physics in the continuum limit, while a 
different phase structure may appear at strong coupling. 
We will not investigate this option here. 

III. SIMULATION ALGORITHM 

Simulations for this work are performed using the 
RHMC algorithm [54] . This is an exact algorithm which 
is flexible enough to cover all the cases of interest with 



good performances, while still maintaining a relatively 
simple form. The algorithm itself is a variation of the 
HMC algorithm [5^ where rational approximations are 
used to compute fractional powers of the fermionic ma- 
trix appearing for a generic number of flavors and 
pseudofermions. Both the HMC and RHMC set up a 
Markov chain of gauge fleld conflgurations yielding the 
desired limiting distribution using the same steps. The 
outline of the algorithms can be summarised as follows: 

i) generate new pseudofermion fields from a heat-bath; 

ii) evolve the gauge field configuration following the flow 
of a fictitious Hamiltonian with randomly chosen initial 
momenta for each link variable - this step is usually 
called molecular dynamic (MD) evolution; iii) perform 
a Metropolis test at the of the MD trajectory to correct 
for errors in the integration of the equation of motion at 
the previous step. 

Our implementation of the HMC /RHMC is a straight- 
forward adaptation of the standard technique to the case 
under consideration. As usual the bottleneck of the al- 
gorithm is the inversion of the Dirac operator that is 
needed during the MD evolution. To reduce the number 
of steps in the integration of the equation of motion and 
thus the number of inversions, we use the second order 
Omelyan integrator [56.] for the MD evolution with differ- 
ent time steps for the gauge and fermion action. Multiple 
pseudofermions and the even-odd preconditioning of the 
fermionic determinant are also used in this paper, but no 
other acceleration techniques are necessary for the pur- 
pose of this work. More extensive simulations will require 
a more careful study of these issues. 

In the rest of this section we will describe the mod- 
ifications required to handle fermions in a generic rep- 
resentation and the implementation details used in our 
code. 



A. HMC molecular dynamics 

In the following the expression for the force needed 
in the Molecular Dynamics is generalized to the case of 
fermions in a generic representation. Since simulations 
in arbitrary representations are still at their early stages, 
we decide to derive the forces in detail, in order to define 
the quantities that appear in the simulations as clearly as 
possible. Introducing for each link variable a conjugate 
momentum in the algebra of the gauge group: 

n{x,^l)^nT%x,^i)T^, (14) 

a fictitious Hamiltonian H is written on the group man- 
ifold as: 

H^Hrr+Hg+Hf, (15) 
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where, assuming Uf = 2: 



X.fl 



(17) 



(18) 



Since the momenta are conjugate to the hnk variables, 
they are defined in the fundamental representation. 

The cases where nf ^ 2 are dealt using the RHMC al- 
gorithm described in the next section. Also note that an 
arbitrary shift s has been included in the fermionic ac- 
tion which will be useful for the discussion of the RHMC 
algorithm. 

Denoting by r the fictitious time of this Hamiltonian 
system, the equation of motion are given by: 



dr 

d 

d^ 



U{x,ii) — 7r(a;, /i)[/(a;, /i) 



(19) 
(20) 



where the RHS of the second equation above defines the 
force F{x,fi) that drives the time evolution of the mo- 
menta. 

For a generic infinitesimal variation of a link variable: 



SU {x, /i) — 5uj{x, /i)[/(x, /x). 



(21) 



where 5uj{x, fi) = iSuj"'{x, n)Tp is an element of the alge- 
bra, F{x, ^) is obtained from the corresponding variation 
of the action through the equations: 

F{x,fi) = Fg{x,fi) + Ff{x,fj.), 
SSg = -iSio,Fg), 
SSf - -{5u;,Ff). 

Since we use link variables in the fundamental represen- 
tation for the gauge action S'g, the variation of the latter 
is the usual one that appears in the HMC molecular dy- 
namics evolution: 

^ -J^T. SLo''ix,^i)Retv[iT^U{x,^i)VHx,^J:^2} 

where V{x, fi) is the sum of the forward and backward 
staples around the link U{x,ij). 

The fermionic force is obtained from the variation of 
the fermionic action as follows. Starting from: 

5Sf = - 0t(Q^ _ s)-i<5(g2„)(Q^ - s)-V, (23) 

let us define: 



r? = {Ql-s)-\^ 

£, = Qmv; 



(24) 
(25) 



using the fact that the matrix (Q^ — s) is hermitean, we 
can rewrite Eq. ([^5]) as 

5Sf = -2Re [C^5{Qm)r]] . (26) 

Inserting the explicit form of Qm, Eq. ([9|), into Eq. (|26|) 
we obtain 

SSf^Re J2 [a^)^SU''{x,fi)^,il^j^)iix + fi)+ 

+77(x)t<5C/«(a:, m)75(1 - 7m)^(^ + A^)] (27) 

where we are implicitly summing over spin and color in- 
dices. We can now write the variation of as: 

SU^ix, fi) = 5uj"-{x, n)U"{x, ^) = i Suj^ix, ti)T^U^{x, n), 

(28) 

where the 6uj°'{x,^) in the above equation are the same 
functions that define the variation of the gauge link in 
Eq. (Ell)- Eq. dlZl) and ^ yield: 



SSf = Suj''ix,fi)Re tre,s [iT^U'^ix, fi)-/5 x 



(29) 

Here the symbol tr^^ indicates the trace over color and 
spin. 

For the sake of convenience we introduce the following 
projectors over the algebra in the representation R: 

Pm = -7lr^<^t^c[tnF] , (30) 

and the following trace operator: 

Tr,^^(77,e) =tr, [j5{l - 1^) {vix + f^) ® ^{x)^ + 

+ax + ^l)(^vix)^}] . (31) 

The forces are then given by: 

Foi^^f^) = -^Pf {U{x,p.)V\x,iJi)) , (32) 
F^{x,^Ji) - (C/«(a;,M)Tr,,^(r;,e)) . (33) 

Note that when R is chosen to be the fundamental 
representation, the usual expressions for the fermionic 
force are recovered. Explicit expressions for the forces 
had already been computed for the case of fermions in 
the adjoint representation of the gauge group SU(2) in 
Ref. [53 • As an analytic check, we have verified that our 
result above agrees with Ec^. (16) in Ref. [svj . 



B. RHMC implementation 

The fermionic part of the HMC hamiltonian, for n f de- 
generate fermions and Npf pseudofermions, can be writ- 
ten as: 



N. 



pf 



(34) 



fc=l 



fc=l 
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For the sake of simplicity we will set all the Ik to be equal: 



Vfc, 



2iV, 



pf 



(35) 



The above decomposition is used in the RHMC algo- 
rithm [s^ l- where rational approximations are used to 
compute the fractional powers of the positive definite 
fermion matrix (5^„. Even though we will work at fixed 
Uf = 2 in this paper, the RHMC is particularly useful if 
one wants an algorithm that can be easily generalized to 
an arbitrary number of fermions. Three different rational 
approximations are used for this implementation. 

The first rational approximation is required in the 
heat-bath update of pseudofermions. In order to gen- 
erate pseudofermions distributed as in Eq. a simple 
two-step process is used. For each pseudofermion we first 
generate a gaussian distributed field (jfk with probability: 



P{(t>k) oc exp[- 



and then we set: 



4>k iQ^n) ' 4>k 



(36) 



(37) 



making use of the fact that (Q^) hermitean (notice the 
plus sign in the exponent.) The RHMC algorithm uses a 
rational approximation to compute the above quantities 
(we need only one approximation since all Ik are equal) : 



di 

E 



The second rational approximation is used to approx- 
imate Eq. p4p during the molecular dynamics evolution 
(as before only one approximation is needed because all 
Ik are equal): 



J2'l'i'^b{Ql)(bk, (39) 

k=l 



Using the formulas derived in Section IHI Al it is easy 
to write the force corresponding to Eq. (|39p . In fact, 
Eq. (j39[) is nothing but a sum of terms of the form 
Eq. (dll) once we put / = 1, s = s^. The RHMC force 
will be then given by a sum over n = 1, . . . , (i2 of terms 
given by Eq. ((33|) multiplied by a factor a\^. It is possible 
to compute all the 77 's defined in Eq. (|24|) corresponding 
to different simultaneously with a multi-shift inverter. 

The third rational approximation is used for the 
Metropolis test. Starting from Eq. (|34]). for each pseud- 
ofermion we can rewrite: 



(41) 



where we used the property that is hermitean. 
rational approximation needed in this case is: 



The 



=1 



the coefficients for the two ap- 
can each be obtained from the 



Notice that if di = 
proximations Tq and rc 
other. 

The coefficients a„, s„ appearing in the rational ap- 
proximations are computed using the Remez algorithm. 
In this implementation we do not compute the coeffi- 
cients "on the fly" but instead the required values are 
taken from a look-up table that has been precomputed, 
according to the following prescription. 

First note that we need to compute rational approxi- 
mations for a function f{x) of the form f{x) — and the 
approximation must be accurate over the spectral range 
of the operator To simplify the computation of the 
table we note that the following proposition holds: if f{x) 
is a homogeneous function of degree I and r{x) is an op- 
timal (in the sense of relative error) rational approxima- 
tion to f{x) over the interval [e,h] to a given accuracy 
then r{kx)/k} is an optimal rational approximation for 
the same function and the same accuracy over the inter- 
val [e/k, h/fc]. Moreover the coefficients of this "rescaled" 
rational approximation are easily obtained from that of 
the original approximation. A simple corollary is that 
we can divide the optimal rational approximations of a 
given homogeneous function f{x) with the same relative 
precision, in classes labeled by the ratio e/h. Within 
each of these classes the coefficients of the rational ap- 
proximations are easily related to each other, so that 
we only need to compute one rational approximation in 
each class. Thus all that is needed in practice is a ta- 
ble containing the coefficients of rational approximations 
for each of these classes, for each function f{x) which 
we want to approximate and for each relative precision 
which is required. At run-time this table is used to gen- 
erate optimal rational approximations rescaling the pre- 
computed coefficients to the desired interval containing 
the spectrum of the matrix Qm • This interval is obtained 
by computing the maximum and minimum eigenvalue of 
on each configuration when needed. 



C. Even-Odd preconditioning 

As already pointed out above, the time required for the 
inversions of the Dirac operator is one of the dominant 
contributions to the total cost of the simulation. The 
convergence of such inversions can be improved using an 
appropriate preconditioning. The idea behind precondi- 
tioning is to rewrite the fermionic determinant as a de- 
terminant (or product of determinants) of one (or more) 
better conditioned matrix (matrices) than the original 
Dirac operator. Very effective preconditionings are at the 
heart of the recent progress in simulations of QCD with 
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dynamical fermions. For the scope of this work, we use 
a simple even-odd preconditioning [s^. We review the 
main features of even-odd preconditioning here in order 
to explain the modifications that are required when con- 
sidering higher-dimensional representations. We divide 
the lattice in a sublattice of even points Ae and another 
sublattice of odd points Aq, and we rewrite the Dirac 
operator Dm as a block matrix: 

™" [Doe 4 + ' ^^^^ 

where each block has a dimension half that of the orig- 
inal Dirac matrix. The diagonal blocks connecting sites 
with the same parity are proportional to the identity ma- 
trix, while off-diagonal blocks connect sites with opposite 
parity. We have (since Dm is 75-hermitean): 

j5Deol5 = Dl . (44) 

The determinant of the Dirac matrix Dm can be rewrit- 
ten as: 

det Dm = det [(4 + - D^oDoe] = det D'^ , (45) 

using the well known expression for the determinant of a 
block matrix. Since the determinant of Dm and of Df" 
are the same the latter can be used in numerical simu- 
lations. Note that the even-odd preconditioned matrix 
only connects sites with the same parity thus it has only 
half of the size of the original Dirac matrix and, like Dm, 
it is 75-hermitean. We define as before the hermitean 
matrix Qm = 75-0^, which will be used in practice. 

The formulation of the HMC algorithm does not 
change and the only difference is that pseudofermions 
fields are now only defined on half of the lattice sites, 
conventionally the even sites in what follows. We now 
derive the explicit expression for the fermionic force for 
the preconditioned system described by the hamiltonian: 

Hf = [(Qm)'-s]"Ve , (46) 

where as before we are assuming Uf = 2 or a rational 
approximation of the actual fractional power function; 
the suffix in is an explicit reminder that the pseud- 
ofermion field is only defined on even sites. Eq. (|26p is 
unchanged : 

SSf = -2Re [?|5(g^°)77e] , (47) 
where as before we have defined: 

Ve = {{Q'^f-sr'clye, (48) 

= QZn.- (49) 

The explicit form of must be used at this point. We 
have: 

= -l5i5DeoDoe + DeoSDoe) . (50) 



Defining 

r]o = DoeVe , (51) 
Co = DoeCe, (52) 

and inserting Eq. (|50p into Eq. (|47l) . we arrive at the 
same expression as before for the variation of the action, 
but with a minus sign: 

6Sf= - 5ZRetr45t/^(x,Ai)Tr,,^(?7,e)] (53) 

From Eq. ([53| and proceeding as before we arrive at the 
final expression for the force. This again coincides with 
Eq. ([55)1 but with the opposite sign: 

F^{x, m) = -P^ {U^'ix, /i)Tr,,^(77, 0) (54) 



D. Some details on the implementation 

As the algorithm is implemented for the first time, we 
describe in this section some general features of our code 
and the technical solutions adopted to manage fermions 
fields defined in an arbitrary representation. 

At this early stage in our physics studies we prioritize 
the flexibility of the code over the optimization aspects. 
The code is written in ANSI C for maximum portability 
to different architectures. This choice is justified, in our 
opinion, by the fact that, on this relatively small scale, 
C codes maintain a much higher overall simplicity and 
clarity compared to e.g. CH — h codes with a comparable 
low overhead. 

The major difference with respect to conventional 
codes - i.e. codes for simulations of the SU(3) gauge 
group with fundamental fermions - is obviously that all 
operations involving gauge links or pseudofermions fields 
must be able to handle matrices and vectors of arbi- 
trary dimension. In the generic case of four-dimensional 
SU(7V) with fermions in a representation R, link matrices 
for the gauge action have dimension N x N, the repre- 
sented gauge links have dimension d^j x df{ and fermion 
fields have dimension Adn.. These numbers indicate how 
the cost of the simulation scales both in memory and 
in the number of floating point operations per applica- 
tion of the Dirac operator as N and R are varied. One 
simple way to handle this objects would be to have func- 
tions that implement operations on them taking as an 
additional parameter the dimensions of the objects them- 
selves. However we found it more flexible to use the fol- 
lowing approach. We use preprocessor macros instead 
of functions^. We fix the gauge group and the represen- 
tation R at compile time as this gives the opportunity 



Inline functions could work just as well, however the handling of 
inlining is a feature which depends on the compiler. 
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to more optimizations (e.g. partial rmrolling of loops)^. 
All of the required macros are automatically generated 
in a pre-compiling step. We use a simple Perl program 
for this purpose whose input parameters are the number 
of colors N and the representation R and whose out- 
put is a C header file containing the macro functions. 
The whole process of header generation and compilation 
is managed through a custom Makefile system and it is 
thus completely transparent to the end user. 

A second difference with respect to conventional codes, 
is that it is convenient, although not strictly necessary, 
to have a second copy of the gauge field in the represen- 
tation R. This is only used in the computation of the 
fermionic force, see Eq. (p3)) . and it is also useful to im- 
pose boundary conditions on the fermion fields. Every 
time the gauge field is updated during the molecular dy- 
namics evolution, the represented gauge field needs to be 
recomputed. 

At this early stage of simulations, we chose to use dou- 
ble precision everywhere in the code to reduce the risk of 
problems arising from numerical accuracy. This is done 
at the expense of optimal performance, which we feel is 
justified in this study of the algorithm implementation 
and behavior. As the source of random numbers, we use 
the RANLUX generator 



IV. BEHAVIOR OF THE ALGORITHM 

In this section we shall discuss some simple tests of the 
RHMC algorithm with the aim to show the correctness 
of our implementation. The majority of the tests per- 
formed for this purpose were run on the lattice T2-B11 
(see Tab. IIV CI below) which corresponds to our biggest 
volume and the lightest mass at /? = 2. For all the tests in 
this section we have used even-odd preconditioning and 
two pseudofermion fields. 
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FIG. 1: Dependence of {exp(Aff)) and (P) on the time-step 
used for the MD integration. 
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A. Preliminary Tests 



FIG. 2: The expectation value (AH) is proportional to (Ar)* 
(upper panel) consistently with the use of a second order in- 
tegrator (the red curve shows the one parameter fit). The ac- 
ceptance probability Pace as measured in the test runs (lower 
panel) is correctly described by the expected large volume 
behavior Pace ~ erfc(-^ {AH)/2) (solid curve, not a fit). 



A simple test of the algorithm is to look at the expec- 
tation value (exp(— Ai?)), where AH is the difference of 
the values of the Hamiltonian H at the beginning and 
at the end of a trajectory in the Molecular Dynamics 
evolution. This expectation value should always be I (a 
result known as Creutz equality [lOl). Moreover since the 
RHMC is an exact algorithm, i.e. there are no corrections 
due to the errors occurring in the Molecular Dynamics in- 
tegration, a good test of the algorithm is to check whether 



^ This also means that recompiling is needed whenever the num- 
ber of colors or the fermions representation is changed. However 
since the total compilation time, including the automatic genera- 
tions of required header, is less than one minute and independent 
on N, this can hardly be a problem. 



or not there is any dependence of any quantity on the in- 
tegration step size At. The average plaquette (P) is a 
good candidate to this purpose because it can be mea- 
sured with high accuracy. 

We show in Fig. [T] the quantities (exp(— Ai/)) (upper 
panel) and the average plaquette (lower panel) for four 
difi^erent step sizes of the MD integration. No statis- 
tically significant deviations from the expected behavior 
are seen. The value of (P) is consistently independent on 
At and the Creutz's equality {exp{~ AH)) = 1 is satis- 
fied. Notice also that for the latter quantity the error bars 
become smaller as At goes to zero, as expected from the 
fact that (AH) is also vanishing in the same limit. This 
is explicitly shown in the upper panel of Fig. (3] As ex- 
pected for a second order integrator (see [6l|), the quan- 
tity {AH) is proportional to At"* with very good accu- 
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FIG. 3: Reversibility test for several values of the time-step 
used for the MD integration. 
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racy (the red curve in the plot is the best fit to that func- 
tional form with one free parameter). The acceptance 
probability Pace can also be used a test of the correct- 
ness of the algorithm. In fact in the large volume limit 
this probability is given by Pace = erfc(^ {^H) /2) [6^ . 
We show the measured value of the acceptance rate as 
a function of the average Hamiltonian variation in the 
lower panel of Fig. [21 The same figure shows also the 
predicted behavior (solid curve). We found a convincing 
agreement with the expectations. 

As a last test in this section we measured reversibility 
violations in the MD. We remind the reader that this is a 
necessary condition for the correctness of the algorithm. 
The discussion will follow that of Ref . [6^ . At the "mi- 
croscopic" level, reversibility violations can occur when 
updating the gauge variables and in the momenta update 
during the MD evolution. We deal with the first source 
of non-reversibility using the trick suggested in Ref. f^: 
the link update Eq. (|19l) is implemented by left multipli- 
cation of an exactly unitary matrix E['k{x, /i)] such that 
E[i:{x^ ^)]E[—'n(x, n)] = 1 and which is also an approxi- 
mation to the exponential map. In this way no possible 
reversibility violations can arise (within machine preci- 
sion) in the link update. For the momenta update, local 
reversibility is guaranteed by the use of double precision 
and the requirement of a small relative residue (10~^ in 
our simulations) for the force calculation. The solution 
of the linear system required for the force calculation 
itself is always started from the same trial solution. At 
the "global" level we measured the reversibility violations 
through the quantity \5H\ defined as the average differ- 
ence between the Hamiltonian of the starting configura- 
tion and the one obtained evolving the system forward 
for a unit of Monte Carlo time and then back (hipping 
the momenta at r = 1) to the original position in phase 
space. This measure of reversibility violations is shown 
in Fig. 131 The global violations to reversibility appear to 
be very small and independent on Ar and hence also on 
(Aif). 



FIG. 4: Time-history for the thermalization of the plaquette 
in the SU(3) theory with two flavours in the fundamental 
representation on a 16* lattice a.t jS = 5.6 and k — 0.15750. 
The black line shows the evolution of the plaquette value using 
our HMC algorithm, while the red line represents the same 
quantity using the DD-HMC algorithm. Trajectory number 
has been scaled by a factor of 7 for the DD-HMC data. 



automatically by a Perl program. As a first non-trivial 
test of our implementation, our code should correctly re- 
produce the results for simulations in the fundamental 
representation of SU(3). The code has been benchmarked 
for SU(3) with Uf — 2 against simulations obtained using 
the DD-HMC algorithm [3] on smaU lattices. Fig. HVBl 
shows the thermalization of the plaquette at (3 = 6.0. 
On small lattices, the DD-HMC only updates a small 
fraction of the links in the system, and therefore the 
thermalization is much slower in units of MC trajecto- 
ries. To compare the two time-histories we have therefore 
rescaled the trajectory number for the DD-HMC data. 
The figure shows a very good agreement between the two 
simulations, in particular for the equilibrium value of the 
plaquette. 

Another test of the algorithm is obtained by simulating 
the SU(3) theory with two flavors in the fundamental 
representation and in the two-index antisymmetric one. 
For the theory with three colors the two representations 
are equivalent and related by charge conjugation. We 
have compared the outcomes of two simulations on 4"^ x 8 
lattice at (3 = 5.6 and k = 0.15600. The time-history 
of the plaquette and its probability distribution for the 
two simulations are summarized in Fig. [5l Once again 
we find a very convincing agreement between the two 
simulations. 



C. SU(2) Uf = 2 two— index symmetric 



B. Checks for SU(3) with n/ = 2 in the 
fundamental representation 

The routines that perform linear algebra operations in 
color space for arbitrary representations are generated 



The theory with two colors and two Dirac fermions in 
the two-index symmetric (25) representation has been 
proposed as a phenomenologically relevant candidate for 
a minimal walking technicolor theory. For the specific 
case of the SU(2) color group, the 25" representation is 
equivalent to the adjoint representation and is therefore 



9 



0.20 r 

0.15 
0.10 
0.05 ■ 
0.00 



SU(3) 2AS 
SU(3) fund 



).55 



0.56 



0.57 0.58 
Plaquette 



0.59 



FIG. 5: Probability distribution of the plaquette. The black 
(respectively red) curve refers to data from a simulation of 
the SU(3) gauge theory on a 4^ x 8 lattice at /3 = 5.6, k — 
0.15600 with two fermions in the fundamental representation 
(respectively the two-index antisymmetric). 
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9.0(9.6) 
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0.6168(73) 


2.6(1.8) 


T2-B11 


16 X 8^ 


0.18797 


1.34 


3840 


0.6347(58) 


13.6(11.5) 



TABLE I: List of runs for the T2 theory at inverse coupling 
/? = 2.0. The standard runs are performed using link vari- 
ables in the real adjoint representation. The primed run is 
performed with link variables in the 25 representation that 
are constructed in an algorithmic way as complex 3x3 ma- 
trices. In order to simplify the comparison with the runs in 
Ref. [231 we list both the values of k and m 



real, which makes it less intensive from the computational 
point of view. Analytical results for the A parameter, 
the renormalization constant, and the possibility of an 
Aoki phase at finite lattice spacing have been presented 
in Ref. [iJ], where the theory was denoted T2. We will 
adopt the same conventions here. Preliminary results 
from lattice simulations of this theory were presented in 
Ref. [i^ . In this subsection we build upon this previous 
experience and provide some more extensive description 
of the behaviour of numerical simulations in the space 
of bare parameters. The list of runs is summarized in 
Tab. IVSlfor inverse coupling l3 = 2.0 and in Tab. HVT^ 
for inverse coupling /? = 2.25. 

This theory allows us to perform a non-trivial check 
of our code generation. The link variables in the higher 
representation can be built in two different manners: in 
one case we construct real matrices according to Eq. ([T]) , 
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FIG. 6: Normalized histogram of the plaquette distribution. 
Data for the SU(2) theory with two flavours in the adjoint 
representation (red curve) are compared with data for two 
flavours in the two-index symmetric representation (black 
curve) . 



while in the other case we use our generic algorithm for 
constructing the 25* complex representation of SU(A^), as 
described in App. |B] As expected, for a given configu- 
ration in the fundamental representation the two repre- 
sentations yield exactly the same matrices; the generic 
algorithm produces complex 3x3 matrices with vanish- 
ing imaginary parts. Evolving the configuration using 
the HMC with the two different representations yields 
compatible results as one can see by the comparison of 
the runs T2-A1 and T2-A1' shown in Fig. WC\ for the 
case of the plaquette. 

The critical value of the hopping parameter has been 
computed in perturbation theory in Ref. [g^l, to which 
we refer the reader for the notation and results used in 
what follows. Using the cactus resummed formula for the 
theory under consideration, we find: 



rUc 



gl2C2{R) X (-0.162857)/co 



(55) 



where cq is a function of defined by solving the follow- 
ing equation: 



Co 



-9o/(16£o) 



1 



245o 



(56) 



For /3 = 2.0 and /3 — 2.25 the above formula yields re- 
spectively rUc = —1.73 [kc = 0.220) and TOc = —1.47 
(kc = 0.198). A non-perturbative determination of the 
chiral limit in the space of bare parameters (/3, k) can be 
obtained by linear extrapolation of the PCAC mass, or of 
the mass squared of the pseudoscalar Goldstone boson. 
Due to the power divergencies in the renormalization of 
the bare quark mass, any perturbative computation of 
TOc is expected to receive non-perturbative corrections 
(9(g-i/9o('')/a). Therefore the perturbative computation 
can at best yield an indication of the location of the crit- 
ical point. The detailed study of these quantities is pre- 
sented in the next section. 

In order to check that the range of k values does not 
lead to exceptionally small values of the eigenvalues of 
the Dirac operator, we have monitored the distribution 
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3 
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0.673737(46) 
0.675184(59) 
0.676649(52) 
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0.2227(10) 
0.07036(90) 
0.05167(50) 
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0.02474(28) 



3.05(36) 
5.9(1.5) 
6.1(1.1) 
4.66(75) 
7.9(1.5) 



0.04436(51) 
0.02836(59) 
0.01520(39) 



3.5(1.5) 
4.2(2.5) 
5.7(3.6) 



TABLE II: List of runs for the T2 theory at inverse coupling /? = 2.25. In the last two columns we give the average of the 
smallest eigenvalue of \Qm \ and its integrated autocorrelation time. 



of the smallest eigenvalue of QJ^ follovifing the study in 
Ref. 3 • The histograms that describe the distribution of 
the lowest eigenvalue of \Qf°\ are displayed in Fig. IIV C] 
for the simulations at (3 — 2.25. At the values of k that 
we have considered, the spectrum of the Dirac operator 
has a clear gap. This is reflected in a smooth behaviour 
of the simulations even at the lightest masses, as can be 
seen in the time-histories for the plaquette and the solver 
number reported in Fig. IIV CI The results obtained with 
fundamental fermions in Ref. Q suggest that the width 
of the distribution scales like a/VV as the continuum and 
thermodynamic limits are approached. A comparison of 
the standard deviations of the eigenvalue distributions for 
the simulations at /? = 2.25 shows that this scaling is well 
verified for the theory under consideration (Fig. IIV Cl) . 

As a final test of the algorithm we monitored the MC 
integration forces for the RHMC algorithm. Compatibly 
with the absence of exceptionally small eignvalues, these 
forces appear to behave smoothly as a function of the 



bare mass and we observe a very mild (consistent with 
zero) dependence on the size of the lattice at constant 
airiQ. In the upper panel of Fig. IIV D I we show the mod- 
ulus of the force |F|.(x,/i)| averaged over all the lattice 
points X and all directions /i for each term of the rational 
approximation indexed by n. Data for two different vol- 
umes at the same value of the bare quark mass are shown 
corresponding to the lattices T2-E2 and T2-F0. The av- 
erage fermionic force appears to change very little once 
the equilibrium is reached which is reflected in the plot 
by the fact that errors on this quantity is really small and 
not visible on the scale of the graph. The realitive error 
for the average force ranges from a few permille for n = 1 
to less than one percent for n = 11. The variation with 
volume is also very small: the two average values appear 
to be always consistent with each other within errors. In 
the bottom panel of Fig. IIV D] we report for completeness 
the size of the rational approximation coefficients used in 
the simulation. Their variation is also quite small, given 
that the distribution of the smallest and largest eigen- 
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FIG. 7: Probability distribution of the lowest eigenvalue of 
the Dirac operator for (3 = 2.25. 




FIG. 8: Time-histories for the plaquette and the solver num- 
ber for the T2-B11 lattice at ^ = 2. 



value of the Dirac operator is quite narrow as discussed 
above. 

The dependence of the average fermionic forces on the 
bare quark mass is shown in Fig. IIVDI for the T2-F lat- 
tice. The change with the quark mass is quite small but 
visible in this case, and more pronounced for higher val- 
ues of n, corresponding to smaller shifts. 



D. Autocorrelation times 

The integrated autocorrelation time for the plaquette 
during our runs is computed in units of molecular dynam- 
ics time following Ref. The results are reported in 
Tab. IIV C[ where an increase of the autocorrelation time 
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FIG. 9: The distribution width AAq of the lowest eigenvalue 
of the Dirac operator scales like V^^^^. Results for /3 — 2.25. 



is clearly visible as we move to lighter fcrmion masses. 
In all cases the autocorrelation time is short compared 
to the total length of the simulations, which guarantees 
that the binned configurations used for our analyses can 
be considered as independent and therefore the statisti- 
cal error is correctly estimated from the variance of the 
ensembles. It is worthwhile to emphasize that extensive 
simulations of QCD have shown that the autocorrelation 
time depends on the trajectory length (66| . and on the 
lattice spacing. Our runs are all performed with trajec- 
tories of unit length. Moreover a critical slowing down of 
the topological modes has been highlighted in the time- 
histories of the topological charge [67] . Therefore we ex- 
pect that the observables that are more strongly corre- 
lated with topology, show a greater correlation length 
as we approach the continuum limit. These issues need 
to be kept in mind as the simulations evolve and more 
precise estimates of the phenomenological quantities will 
be sought. Finally, one should keep in mind that auto- 
correlations depend on the observables, and should be 
monitored for all the relevant quantities in order to fully 
control the systematics in the simulations. For our runs 
at /3 = 2.25 we changed the MC integrator parameters 
to keep the plaquette autocorrelation time roughly con- 
stant, which is why the data in Tab. IIV CI do not show 
the same increase in autocorrelation times. 



V. BENCHMARK MESONIC OBSERVABLES 

Phenomenological results on the spectrum of the T2 
theory are extracted from the study of lattice two-point 
functions. In this study we do not attempt to perform an 
investigation of the phenomenology of this theory, which 
will be presented in a forthcoming paper. Our aim is to 
provide lattice results useful as a benchmark for future 
investigations, and to define the observables and the anal- 
ysis procedure that we use here and in subsequent papers. 
In this section we describe in detail the procedure used 
to extract mesonic observables, such as the mass of the 
pseuodoscalar or vector mesons, and the corresponding 
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FIG. 10: Volume dependence of the average MC forces for 
each term in the rational approximation (upper panel) and 
average value of the approximation coefficients (lower panel). 
Data are obtained on the T2-E2 and T2-F0 lattices. Error 
bars in all cases are much smaller than the size of the symbols. 
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FIG. 11: Mass dependence of the average MC fermionic forces 
for each term in the rational approximation. Data are ob- 
tained on the T2-F lattices. Error bars are not shown as they 
are much smaller than the symbols used in the plot. 



results on small lattices at (3 — 2.0. 



A. Definitions 

Let r and V be two generic matrices in the Clifford 
algebra, we define the two-point correlator at zero mo- 
mentum as follows: 

frr'it) = ^(^i(x,t)r^2(x,t)^2(0)rVi(0)) , (57) 



where "01 and "02 represent two different flavors of de- 
generate fermion fields, so that we only consider flavor 



non-singlet bilinears. Denoting the space-time position 
(x, t) by X and performing the Wick contractions yields: 

{Mx)TMx)4'2m'M0)) = -tT[rs{x)r's{~x)] 

= -tr [TS{x)r'^5S\xh5] 

In practice we invert the Hermitean Dirac operator Q = 
75D by solving the equation: 



Q{x,y)ABVB'°iy) 



^A.A^x,0 



(58) 



where capital Latin letters like A = {a, a} are collective 
indices for color and spin, and A, x = is the position of 
the source for the inverter. The inversion is performed us- 
ing a QMR recursive algorithm with even-odd precondi- 
tioning of the Dirac operator, which is stopped when the 
residue is less than 10~*. Using the solution ry obtained 
from the inversion, the above correlator is re-expressed 
as: 



(...) 



D.y/ \>t 

cdVa 



(59) 



where f = 75r, andj" = 75r'. 
Following Ref. 



68| 



masses and decay constants for 
the pseudoscalar meson are extracted from the asymp- 
totic behaviour of the correlators /pp and /ap at large 
Euclidean time. The pseudoscalar mass and the vacuum- 
to-meson matrix element are obtained from the correla- 
tor of two pseudoscalar densities: 



fpp{t) = 



f^2 



PS 



exp [-Mpst] + 



(60) 



The meson mass is obtained by fitting the effective mass 
to a constant, while the coupling Gps is extracted from 
the amplitude of the two-point function fpp. For the 
precise definition of the effective mass and coupling, we 
refer the reader to Ref. [63 • On the other hand the ratio 



1 



(61) 



yields the PC AC mass m with corrections of 0{a) for 
the unimproved theory. On the relatively small lattices 
that we have used in this study, it is difficult to isolate 
clearly the contribution from the lowest state, that dom- 
inates the large-time behaviour of two-point correlators. 
Fig. [12] illustrates the typical quality of the plateau that 
is fitted to extract the pion mass. The data on the plot 
correspond to the effective mass obtained from two fpp 
two-point function, whose asymptotic behaviour is gov- 
erned by the mass of the lightest pseudoscalar meson. 
Lattices with a smaller time extent do not show such a 
clear plateau; for these smaller lattices the determination 
of the masses is affected by larger systematics due to the 
contamination from heavier states. 

Note that the decay constant is not computed directly; 
it is obtained from the values computed above as: 



-Fps 



Gps- 



(62) 



13 



0.36 



0.33 




lattice 



aMps 



a^Gps 



aMy My /Ma 



FIG. 12: Effective mass plot for the pion mass. The points 
correspond to the effective mass extracted from fpp{t). The 
data refer to the T2-F0 lattice. 



lattice 



aM-ps 



a^Gps 



aFps 



T2-A1 
T2-A2 
T2-A3 
T2-A4 
T2-A5 
T2-A6 
T2-A7 
T2-A8 



1.0199(14) 
0.7997(14) 
0.6718(14) 
0.5188(13) 
0.3129(23) 
0.2599(39) 
0.1601(30) 
0.0921(16) 



2.5711(12) 
2.2313(16) 
2.0212(19) 
1.7708(24) 
1.362(7) 
1.196(18) 
0.868(27) 
0.677(24) 



1.134(8) 
1.277(10) 
1.351(12) 
1.435(14) 
1.356(34) 
1.157(59) 
0.776(47) 
0.607(27) 



2.5814(13) 
2.2516(17) 
2.0517(21) 
1.8222(26) 
1.457(8) 
1.280(20) 
0.901(31) 
0.684(28) 



0.3499(11) 
0.4102(16) 
0.4442(20) 
0.4748(26) 
0.4578(68) 
0.421(12) 
0.3301(97) 
0.2433(61) 



TABLE III: Results for the PCAC mass, for the pseudoscalar 
mass and coupling, for the vector meson mass and the pseu- 
doscalar decay constant, from simulations on the T2-A lat- 
tices. 



The decay constant extracted from bare lattice correla- 
tors is related to its continuum counterpart by the renor- 
malization constant Za which has been recently com- 
puted in perturbation theory in Ref. [64]. 

Finally the mass of the vector state is extracted from 
the /vv correlator, again following the procedure out- 
lined in Ref. M. 



B. Results 

A first set of results were obtained from runs on the 
T2-A lattices. Our results are summarized in Tab. IIIII 
The results obtained on such small lattices are affected 
by large systematic errors, and hence are not suitable for 
reliable phenomenology. On one hand the limited exten- 
sion of the T2-A lattices in the time direction means that 
it is virtually impossible to identify a proper plateau. As 
a consequence we estimate the mass of the relevant states 
from the value of the effective mass at the centre of the 
lattice. For the heavier masses we find that our results 



T2-B7 0.2209(30) 1.149(11) 1.190(56) 1.269(12) 1.51(22) 
T2-B8 0.1874(40) 1.044(18) 1.032(73) 1.163(20) 1.75(26) 
T2-B9 0.1624(37) 0.952(19) 0.875(63) 1.067(21) 1.77(15) 
T2-B10 0.1307(42) 0.838(29) 0.714(66) 0.952(30) 1.50(18) 
T2-B11 0.0455(14) 0.356(34) 0.214(12) 0.366(43) 0.67(25) 

TABLE IV: Fitted values for the masses and vacuum-to- 
meson matrix element in lattice units for the T2-B lattices. 



lattice 



aFp 



My/Fp 



T2-B7 0.2209(30) 5.98(14) 0.399(14) 3.19(12) 

T2-B8 0.1874(40) 5.82(24) 0.354(20) 3.29(19) 

T2-B9 0.1624(37) 5.58(26) 0.314(18) 3.40(20) 

T2-B10 0.1307(42) 5.30(41) 0.266(21) 3.58(30) 

T2-B11 0.0455(14) 2.79(39) 0.154(8) 2.38(31) 

TABLE V: Fitted values for the relevant combinations of m, 
Mps, Gps, and My for the T2-B lattices. The value of the 
PCAC mass is reported again for clarity. 



are in agreement with the data presented in Ref. psj]'^. 
However, for these larger values of the bare mass, we 
see that all masses are 0(1) in units of the inverse lat- 
tice spacing, and therefore we expect these results to be 
affected by large lattice artefacts. Smaller masses are 
needed in order to identify the chiral dynamics that is rel- 
evant for phcnomcnological studies. The lattices T2-A6, 
T2-A7, and T2-A8 yield smaller masses, but the light- 
est pseudoscalar meson has a mass aMps = 0.68, which 
is still large in units of the UV cutoff. The study of the 
eigenvalue distributions presented in the previous section 
suggests that it is not possible to go to lighter masses on 
the T2-A lattices without entering the regime where the 
algorithm becomes unstable, or the system gets close to 
a phase transition. 

The results obtained on the T2-B lattices for the pseu- 
doscalar mass, the PCAC mass, the vector meson mass, 
the axial vector meson mass, and the vacuum-to-meson 
matrix element are reported in Tab. IIVI As a result of 
the reduced width in the eigenvalue distribution, we can 
afford to simulate closer to the chiral limit, with a PCAC 
mass smaller than 0.1. 

Several combinations of interest are also summarized 
in Tab. El the ratio of the pseudoscalar mass squared to 
the PCAC mass, the bare pseudoscalar decay constant, 
and the mass of the vector meson in units of the pseu- 
doscalar decay constant. They are computed from the 
primary observables discussed above, and the error prop- 
agation is done by the usual jackknife method. The main 



^ The reader may notice however that a different normalization 
for aFps has been used in Ref. |25l|. The two choices differ by 
a factor -v/S. Our choice for the normalization of aFpg yields 
F„ = 93 MeV in QCD. 
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phenomenological conclusions of this paper are obtained 
from these values. 

The computation of the pseudoscalar decay constant in 
the chiral limit is the method of choice to set the lattice 
scale. In a technicolor theory the decay constant is re- 
lated to the vacuum expectation value of the Higgs field, 
and therefore Fps — 250 GeV. A realistic determination 
of the physical value of Fps would be beyond the scope 
of this work. Here we simply report the values of Fps in 
Tab. Inland El 

The ratio My/Fpg is also an interesting quantity of a 
Technicolor theory. Our result for this quanity are sum- 
marized in Tab. Ivl 

We stress again that the results in this section are only 
benchmarks: they are useful for checking future numeri- 
cal simulations but not to draw any solid phenomenolog- 
ical conclusion about this theory. 



VI. CONCLUSIONS 

In this paper we have described the generalization of 
the HMC and RHMC algorithms to the case of generic 
representations of the color group SU(A^). We have de- 
veloped a general framework that can deal with Wilson 
fermions in arbitrary representations and generic num- 
ber of colors N. We have put the emphasis in describing 
in detail the algorithmic issues that need to be faced in 
generalizing existing code to arbitrary representations. 
Since simulations for these new theories are still in their 
infancy, and even the simplest results are not known a 
priori, numerous independent studies will be welcome. 
We have provided an extensive number of tests of the 
algorithm and a detailed study of its behavior for the 
SU(2) gauge theory with fermions in the adjoint repre- 
sentation. Our benchmark result are ideally suited for 
validating lattice simulations of theories with fermions in 
higher representations. 

We have presented benchmark results on small 8"^ x 16 
lattices, giving full the details of the procedure we used, 
so as to make all the results presented in this work re- 
producible beyond any ambiguity. On these small lattice, 
our results should be easily reproducible without any ma- 
jor investement of computational resources. We provided 
a large number of benchmark quantities, ranging from the 
average value of the plaquette to mesonic observables. 

The results in this paper are the first step in a 
more comprehensive program that aims to study nonper- 
turbative phenomena beyond QCD. Robust results for 
the spectrum and decay constants in technicolor theo- 
ries are an important ingredient to search for strongly- 
interacting BSM physics at the LHC. The techniques that 
have been developed for QCD offer all the theoretical 
and algorithmic tools needed to develop a comprehen- 
sive study of technicolor on the lattice. Lattice results 
will be important in order to test these theories as po- 
tential candidates for BSM physics. We stress again that 
an uncompromising theoretical formulation is mandatory 



for such studies that are investigating unknown theories 
in order to have some real predictive power. First nu- 
merical results should be benchmarked carefully against 
each other in order to verify that systematics are under 
control. 

It is worthwhile to emphasize once again that the al- 
gorithmic framework that we have developed here can be 
readily used to study the planar orientifold equivalence 
and lattice supersymmetry. Preliminary results in this 
direction have already appeared jl4| . 
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Appendix A: Group— theoretical factors 

The normalization of the generators in a generic rep- 
resentation R of SU(A^) is fixed by requiring that: 

[T^,T},]^ir'^n, (Ai) 

where the structure constants /"'"^ are the same in all 
representations. We define: 

ti-R [t'^T'') = tr (T]^T]^) = (A2) 
E(^fl^fl)AB = C2{R)5ab. (A3) 

a 

and hence: 

TR = j^:^C2{R)dR (A4) 

where dp is the dimension of the representation R. The 
quadratic Casimir operators may be computed from the 
Young tableaux of the representation of SU(Af) by using 
the formula: 

1 / ™ 2 \ 

C2{R)^^{'nN + Y,nr{n^ + l-2i)-'^\ (A5) 
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R 




Tr 


C2iR) 


fund 


iV 


1/2 


(N^ - 1)/{2N) 


Adj 


- 1 


TV 


N 


2S 


N{N + l)/2 


(iV + 2)/2 


C2{F)2(N + 2)/{N + 1) 


2AS 


N{N - l)/2 


(7V-2)/2 


C2{F)2{N -2)/{N -I) 



TABLE VI: Group invariants used in this work 



where n is the number of boxes in the diagram, i ranges 
over the rows of the Young tableau, m is the number of 
rows, and rii is the number of boxes in the i-th row. 

The quantities du, Tr, C2{R) are hsted in Tab. IVll 
for the fundamental, adjoint, 2-index symmetric, and 2- 
indcx antisymmetric representations. 



Appendix B: Two— index representations 

The hermitean generators Tp for the fundamental rep- 
resentation used are of the following form. For each pair 
of integers l<k<l<N,we define two generators as: 



) , (Bi) 

{Tp' )mn ~ -^{SrakSnl ^ Smldnk) , (B2) 



and for each 1 < k < N one more generator is defined 
as: 



1 



A/2fc(fc+ 1 



:diag(l, 1, 



-fc,0, ...,0) . (B3) 



k+l terms 



The generators Tp are normalized so that Tp — 1/2. 
The generators for the other representations are obtained 
as follows. 

We first give the explicit form for the representation 
functions R which map U — ^ U^. We define for each 
representation an orthonormal base br for the appropri- 
ate vector space of matrices. 

For the Adjoint representation we define the base 
BAdj for the N X N traceless hermitean matrices to be 



e^^^. = T^/y/T^, a = 1,...,N'^ - 1 (i.e. proportional 
to the generators of the fundamental representation and 
normalized to 1.) 

For the two-index Symmetric representation the base 
e^'-'\ with i < j, for the N x N symmetric matrices is 
given by: 



i^j, {4'^)kl = ^{S^kSJl + SJkS^l) , (B4) 



1 



/ (fi) N 



V2 



(B5) 



For the two-index Antisymmetric representation the 

(ii) 

base e){g , with i < j, for the N x N antisymmetric ma- 
trices is given by: 

(B6) 



{e-^s^)ki = -^{^ikSji - SjkSii) . 



The maps R are explicitly given by: 

a,6 = l,...,A^2_;L, (B7) 



{R^U) 



'^{i3)ilk) 



tr 



tr 



«<i,^<fc, (B8) 
i<j,l<k. (B9) 



The generators used are defined as the image of 
the generators in the fundamental under the differential 
of the maps R defined above: Tjj = R^Tp. Explicit 
expression can easily be worked out form the definition 
above. The invariants Tr and C2{R) for the generators 
defined in this way are given in Tab. I VII 
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